SaoMNK Tutorial: Exogenous Shocks
S. Downing
2025-02-13
This project begins from the simple premise that the task of
endogenizing the NK search landscape can be done by representing the
search environment as a binary bipartite network of M
actors affiliating with N components. This bipartite
network can be then be analyzed according to the Stochastic
Actor-Oriented Model (SAOM) (Snijders,
1996).
Thus, the Stochastic Actor-Oriented MNK model, abbrivated
SaoMNK, is designed for running strategic search
simulation, testing, and experimentation by leveraging the
RSiena package, an R implementation of
Simulation Investigation for Empirical Network Analysis (SIENA).
This tutorial offers a basic introduction to SaoMNK.
############### Load R6 Class DEPENDENCIES ############################
## Biparite Environment Search Simulation Class
SaomNkRSienaBiEnv <- source(file.path(dir_r, 'SAOM_NK_R6_model.R'))$value
## default settings: Users do not change; TODO: implment within restricted class attributes
DV_NAME <- 'self$bipartite_rsienaDV'
## define simulation iterations based on
steps_per_actor <- 121. Baseline Case: No Shocks
1.1 Environment Configuration
#
environ_params <- list(
M = 12, ## Actors
N = 12, ## Components
BI_PROB = 0, ## Environmental Density (DGP hyperparameter)
component_matrix_start = 'rand', ##**TODO** Implement: 'rand','modular','semi-modular',...
rand_seed = 1234,
plot_init = F,
name = '_test_tutorial_nb_'
)
#
env1 <- SaomNkRSienaBiEnv$new(environ_params)##
## TEST FROM CALLED CLASS INIT *BEFORE* BASE INIT
##
## CALLED _BASE_ INIT
##
## TEST FROM CALLED CLASS INIT *AFTER* BASE INIT1.2 Structrual Model
SAOM Objective Function serves as the stochastic actor’s utility function for strategic search.
#
strategies <- list(
egoX = c(-1, 0, 1),
inPopX = c( 1, 0, -1)
)
## 2.b. Component Payoffs vector
set.seed(12345)
component_payoffs <- runif(environ_params$N, min = 0, max = 1)
## 2. Strategies sets the objective function as a linear combination of network stats across DVs
#
actor_strats_list <- lapply(strategies, function(strat) rep(strat, environ_params$M/length(strat)) )
component_int_mat <- create_block_diag(environ_params$N, round(environ_params$N/4))
# dyad_cov_XWX <- ( outer(component_payoffs, component_payoffs, '*') *
# create_block_diag(environ_params$N, round(environ_params$N/3))
# )
dyad_cov_X <- matrix(runif(n = environ_params$N*environ_params$M, min = -.5, max= 1), nrow=environ_params$M)
# dyad_cov_X <- matrix(runif(environ_params$N*environ_params$M) - runif(environ_params$N*environ_params$M),
# nrow=environ_params$M)
# dyad_cov_XWX_X <- t( dyad_cov_XWX %*% t(dyad_cov_X) )
#
structure_model <- list(
dv_bipartite = list(
name = 'self$bipartite_rsienaDV',
rates = list(
list(effect='Rate', parameter= .111, dv_name=DV_NAME, fix=FALSE )#, ##interaction1 = NULL
),
effects = list( ##**STRUCTURAL EFFECTS -- dyadic/network endogeneity sources**
list(effect='density', parameter= 0, dv_name=DV_NAME, fix=T ), ##interaction1 = NULL
list(effect='outAct', parameter= 0, dv_name=DV_NAME, fix=T ),
list(effect='inPop', parameter= 0, dv_name=DV_NAME, fix=T ) #interaction1 = NUL
),
## COVARIATE EFFECTS
coCovars = list(
##** COMPONENTS : MONADIC CONSTANT COVARIATE EFFECTS **##
list(effect='altX', parameter= 0, dv_name=DV_NAME, fix=T,
interaction1='self$component_1_coCovar', x = component_payoffs
),
# list(effect='outActX', parameter= 0, dv_name=DV_NAME, fix=T,
# interaction1='self$component_1_coCovar', x = component_payoffs
# ),
##** STRATEGIES : MONADIC CONSTANT COVARIATE EFFECTS **##
list(effect='egoX', parameter= .5, dv_name=DV_NAME, fix=T,
interaction1='self$strat_1_coCovar', x = actor_strats_list[[1]]
), #interaction1 = NULL
list(effect='inPopX', parameter= 1, dv_name=DV_NAME, fix=T,
interaction1='self$strat_2_coCovar', x = actor_strats_list[[2]]
)
),
##**MONADIC TIME-VARYING COVARIATE EFFECTS -- DYNAMIC STRATEGY PROGRAMS**
varCovars = list(),
##**DYADIC CONSTANT COVARIATE EFFECTS -- EXOGENOUS INTERACTION MATRIX**
coDyadCovars = list(
list(effect='XWX', parameter= .1, dv_name=DV_NAME, fix=T, nodeSet=c('COMPONENTS','COMPONENTS'),##if M=N, must provide nodeSet
interaction1='self$component_1_coDyadCovar',
x = component_int_mat ## component-[actor]-component dyads
) #,
# list(effect='X', parameter= 0, dv_name=DV_NAME, fix=T, nodeSet=c('ACTORS','COMPONENTS'), ##if M=N, must provide nodeSet
# interaction1='self$component_2_coDyadCovar',
# x = dyad_cov_X ## deltas = changes of payoff contributions from each actor-component
# )
),
##**DYADIC TIME-VARYING COVARIATE EFFECTS -- DYNAMIC INTERACTION MATRIX**
varDyadCovars = list(),
interactions = list(
# list(effect='egoX|XWX', parameter= 0, dv_name=DV_NAME, fix=T,
# interaction1='self$strat_1_coCovar',
# interaction2='self$component_1_coDyadCovar'
# ),
# list(effect='inPopX|altX', parameter= 0, dv_name=DV_NAME, fix=T,
# interaction1='self$strat_2_coCovar',
# interaction2='self$component_1_coCovar'
# ),
# list(effect='egoX|altX', parameter= 0, dv_name=DV_NAME, fix=T,
# interaction1='self$strat_2_coCovar',
# interaction2='self$component_1_coCovar'
# )
)
)
)env1$preview_effects(structure_model, filter=FALSE)## Effects documentation written to file C:/Users/sdr8y/OneDrive - University of Missouri/Research/Search_networks/SaoMNK/R/_rsiena_effects_doc_.html .# ## Uncomment for HTML output with filterable data table
# env1$preview_effects(structure_model, filter=TRUE)1.3 Run RSiena Search Process
## TODO: PICK UP WITH coDydCovar Interation Matrix
## Run Rsiena search using variable parameters in theta_matrix
env1$search_rsiena(
structure_model,
iterations = env1$M * steps_per_actor,
digits = 4,
run_seed = 12345
)##
## Rate effects i=1, j=1
## $effect
## [1] "Rate"
##
## $parameter
## [1] 0.111
##
## $dv_name
## [1] "self$bipartite_rsienaDV"
##
## $fix
## [1] FALSE
##
## effectName include fix test initialValue
## 1 basic rate parameter self$bipartite_rsienaDV TRUE FALSE FALSE 0.1
## parm
## 1 0
## effectName include fix test initialValue
## 1 basic rate parameter self$bipartite_rsienaDV TRUE FALSE FALSE 0
## parm
## 1 0.111
##
## structural effects i=1, j=1
## $effect
## [1] "density"
##
## $parameter
## [1] 0
##
## $dv_name
## [1] "self$bipartite_rsienaDV"
##
## $fix
## [1] TRUE
##
## effectName include fix test initialValue parm
## 1 outdegree (density) TRUE TRUE FALSE -1.60944 0
## effectName include fix test initialValue parm
## 1 outdegree (density) TRUE TRUE FALSE 0 0
##
## structural effects i=1, j=2
## $effect
## [1] "outAct"
##
## $parameter
## [1] 0
##
## $dv_name
## [1] "self$bipartite_rsienaDV"
##
## $fix
## [1] TRUE
##
## effectName include fix test initialValue parm
## 1 outdegree - activity TRUE TRUE FALSE 0 0
## effectName include fix test initialValue parm
## 1 outdegree - activity TRUE TRUE FALSE 0 0
##
## structural effects i=1, j=3
## $effect
## [1] "inPop"
##
## $parameter
## [1] 0
##
## $dv_name
## [1] "self$bipartite_rsienaDV"
##
## $fix
## [1] TRUE
##
## effectName include fix test initialValue parm
## 1 indegree - popularity TRUE TRUE FALSE 0 0
## effectName include fix test initialValue parm
## 1 indegree - popularity TRUE TRUE FALSE 0 0
##
## coCovars i=1, j=1
## $effect
## [1] "altX"
##
## $parameter
## [1] 0
##
## $dv_name
## [1] "self$bipartite_rsienaDV"
##
## $fix
## [1] TRUE
##
## $interaction1
## [1] "self$component_1_coCovar"
##
## $x
## [1] 0.72090390 0.87577319 0.76098233 0.88612457 0.45648096 0.16637179
## [7] 0.32509539 0.50922434 0.72770525 0.98973694 0.03453544 0.15237349
##
## effectName include fix test initialValue parm
## 1 self$component_1_coCovar alter TRUE TRUE FALSE 0 0
## effectName include fix test initialValue parm
## 1 self$component_1_coCovar alter TRUE TRUE FALSE 0 0
##
## coCovars i=1, j=2
## $effect
## [1] "egoX"
##
## $parameter
## [1] 0.5
##
## $dv_name
## [1] "self$bipartite_rsienaDV"
##
## $fix
## [1] TRUE
##
## $interaction1
## [1] "self$strat_1_coCovar"
##
## $x
## [1] -1 0 1 -1 0 1 -1 0 1 -1 0 1
##
## effectName include fix test initialValue parm
## 1 self$strat_1_coCovar ego TRUE TRUE FALSE 0 0
## effectName include fix test initialValue parm
## 1 self$strat_1_coCovar ego TRUE TRUE FALSE 0 0.5
##
## coCovars i=1, j=3
## $effect
## [1] "inPopX"
##
## $parameter
## [1] 1
##
## $dv_name
## [1] "self$bipartite_rsienaDV"
##
## $fix
## [1] TRUE
##
## $interaction1
## [1] "self$strat_2_coCovar"
##
## $x
## [1] 1 0 -1 1 0 -1 1 0 -1 1 0 -1
##
## effectName include fix test initialValue
## 1 ind. pop.^(1/#) weighted self$strat_2_coCovar TRUE TRUE FALSE 0
## parm
## 1 1
## effectName include fix test initialValue
## 1 ind. pop.^(1/#) weighted self$strat_2_coCovar TRUE TRUE FALSE 0
## parm
## 1 1
##
## coDyadCovars i=1, j=1
## $effect
## [1] "XWX"
##
## $parameter
## [1] 0.1
##
## $dv_name
## [1] "self$bipartite_rsienaDV"
##
## $fix
## [1] TRUE
##
## $nodeSet
## [1] "COMPONENTS" "COMPONENTS"
##
## $interaction1
## [1] "self$component_1_coDyadCovar"
##
## $x
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 1 1 1 1 0 0 0 0 0 0 0 0
## [2,] 1 1 1 1 0 0 0 0 0 0 0 0
## [3,] 1 1 1 1 0 0 0 0 0 0 0 0
## [4,] 1 1 1 1 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 1 1 1 1 0 0 0 0
## [6,] 0 0 0 0 1 1 1 1 0 0 0 0
## [7,] 0 0 0 0 1 1 1 1 0 0 0 0
## [8,] 0 0 0 0 1 1 1 1 0 0 0 0
## [9,] 0 0 0 0 0 0 0 0 1 1 1 1
## [10,] 0 0 0 0 0 0 0 0 1 1 1 1
## [11,] 0 0 0 0 0 0 0 0 1 1 1 1
## [12,] 0 0 0 0 0 0 0 0 1 1 1 1
##
## effectName include fix test initialValue
## 1 XW=>X closure of self$component_1_coDyadCovar TRUE TRUE FALSE 0
## parm
## 1 0
## effectName include fix test initialValue
## 1 XW=>X closure of self$component_1_coDyadCovar TRUE TRUE FALSE 0
## parm
## 1 0.1
##
##
## self$rsiena_data :
##
## Dependent variables: self$bipartite_rsienaDV
## Number of observations: 2
##
## Nodesets ACTORS COMPONENTS
## Number of nodes 12 12
##
## Dependent variable self$bipartite_rsienaDV
## Type bipartite
## Observations 2
## First nodeset ACTORS
## Second nodeset COMPONENTS
## Densities NA NA
##
## Constant covariates: self$component_1_coCovar, self$strat_1_coCovar, self$strat_2_coCovar
## Constant dyadic covariates: self$component_1_coDyadCovar
## [1] "density"
## [1] "inPop"
## [1] "outAct"
## [1] "XWX"
## [1] "altX"
## [1] "egoX"
## [1] "inPopX"
##
##
## theta_matrix summary (first and last 5 rows) :
##
## density inPop outAct XWX altX egoX inPopX
## 1 0 0 0 0.1 0 0.5 1
## 2 0 0 0 0.1 0 0.5 1
## 3 0 0 0 0.1 0 0.5 1
## 4 0 0 0 0.1 0 0.5 1
## 5 0 0 0 0.1 0 0.5 1
## 140 0 0 0 0.1 0 0.5 1
## 141 0 0 0 0.1 0 0.5 1
## 142 0 0 0 0.1 0 0.5 1
## 143 0 0 0 0.1 0 0.5 1
## 144 0 0 0 0.1 0 0.5 1
## If you use this algorithm object, siena07 will create/use an output file C:/Users/sdr8y/OneDrive - University of Missouri/Research/Search_networks/SaoMNK/R/_test_tutorial_nb__173944254788.txt .
##
## Start phase 0
## theta: 0 0 0 0 0 0 0
##
## Start phase 3
## Phase 3 Iteration 100 Progress 69%
## Parameter values used for simulations
##
## Mean Standard
## value Deviation
##
## Rate parameters:
## 0 Rate parameter NA ( NA )
##
## Other parameters:
## 1. eval outdegree (density) 0.0 ( NA )
## 2. eval indegree - popularity 0.0 ( NA )
## 3. eval outdegree - activity 0.0 ( NA )
## 4. eval XW=>X closure of self$component_1_coDyadCovar 0.1 ( NA )
## 5. eval self$component_1_coCovar alter 0.0 ( NA )
## 6. eval self$strat_1_coCovar ego 0.5 ( NA )
## 7. eval ind. pop.^(1/1) weighted self$strat_2_coCovar 1.0 ( NA )
##
## Simulated means and standard deviations
## 1. Number of ties 0.882 0.324
## 2. Sum of squared indegrees 0.882 0.324
## 3. Sum of squared outdegrees 0.882 0.324
## 4. XW=>X closure of self$component_1_coDyadCovar 0.000 0.000
## 5. Sum of indegrees x self$component_1_coCovar -0.017 0.307
## 6. Sum of outdegrees x self$strat_1_coCovar -0.076 0.785
## 7. indegree pop.^(1/1) weighted self$strat_2_coCovar 0.076 0.785
##
##
## Simulated statistics are in x$sf
## and simulated dependent variables in x$sims, where x is the created object.
##
## Total of 144 iteration steps.
##
## Covariance matrix of estimates (correlations below diagonal)
##
## 0 0 0 0 0 0 0
## NaN 0 0 0 0 0 0
## NaN NaN 0 0 0 0 0
## NaN NaN NaN 0 0 0 0
## NaN NaN NaN NaN 0 0 0
## NaN NaN NaN NaN NaN 0 0
## NaN NaN NaN NaN NaN NaN 0
##
## Derivative matrix of expected statistics X by parameters:
##
## 0.103 0.103 0.103 0.000 -0.001 -0.031 0.031
## 0.103 0.103 0.103 0.000 -0.001 -0.031 0.031
## 0.103 0.103 0.103 0.000 -0.001 -0.031 0.031
## 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## -0.002 -0.002 -0.002 0.000 0.093 -0.016 0.016
## -0.035 -0.035 -0.035 0.000 -0.001 0.055 -0.055
## 0.000 0.000 0.000 0.000 0.000 0.000 0.000
##
## Covariance matrix of X (correlations below diagonal):
##
## 0.105 0.105 0.105 0.000 -0.002 -0.009 0.009
## 1.000 0.105 0.105 0.000 -0.002 -0.009 0.009
## 1.000 1.000 0.105 0.000 -0.002 -0.009 0.009
## NaN NaN NaN 0.000 0.000 0.000 0.000
## -0.020 -0.020 -0.020 NaN 0.094 -0.016 0.016
## -0.036 -0.036 -0.036 NaN -0.067 0.617 -0.617
## 0.036 0.036 0.036 NaN 0.067 -1.000 0.617
##
##
##
## Simulated Decision Chain Summary:
##
## dv_type dv_type_bin dv_varname id_from
## Length:144 Min. :0 Length:144 Min. : 1.000
## Class :character 1st Qu.:0 Class :character 1st Qu.: 4.000
## Mode :character Median :0 Mode :character Median : 7.000
## Mean :0 Mean : 6.861
## 3rd Qu.:0 3rd Qu.:10.000
## Max. :0 Max. :12.000
## id_to beh_difference reciprocal_rate LogOptionSetProb
## Min. : 1.000 Min. :0 Min. :0.08333 Min. :-2.485
## 1st Qu.: 4.000 1st Qu.:0 1st Qu.:0.08333 1st Qu.:-2.485
## Median : 8.000 Median :0 Median :0.08333 Median :-2.485
## Mean : 7.486 Mean :0 Mean :0.08333 Mean :-2.485
## 3rd Qu.:11.000 3rd Qu.:0 3rd Qu.:0.08333 3rd Qu.:-2.485
## Max. :13.000 Max. :0 Max. :0.08333 Max. :-2.485
## LogChoiceProb diagonal stability tie_change
## Min. :-3.034 Length:144 Mode :logical Mode :logical
## 1st Qu.:-2.614 Class :character FALSE:127 FALSE:17
## Median :-2.565 Mode :character TRUE :17 TRUE :127
## Mean :-2.555
## 3rd Qu.:-2.534
## Max. :-2.114
## chain_step_id
## Min. : 1.00
## 1st Qu.: 36.75
## Median : 72.00
## Mean : 72.37
## 3rd Qu.:108.25
## Max. :144.00
## [1] 144 13
## dv_type dv_type_bin dv_varname id_from id_to beh_difference
## 1 Network 0 self$bipartite_rsienaDV 11 10 0
## 2 Network 0 self$bipartite_rsienaDV 4 7 0
## 3 Network 0 self$bipartite_rsienaDV 12 1 0
## 4 Network 0 self$bipartite_rsienaDV 9 1 0
## 5 Network 0 self$bipartite_rsienaDV 6 5 0
## reciprocal_rate LogOptionSetProb LogChoiceProb diagonal stability tie_change
## 1 0.08333333 -2.484907 -2.564949 FALSE FALSE TRUE
## 2 0.08333333 -2.484907 -2.613646 FALSE FALSE TRUE
## 3 0.08333333 -2.484907 -2.534215 FALSE FALSE TRUE
## 4 0.08333333 -2.484907 -2.534215 FALSE FALSE TRUE
## 5 0.08333333 -2.484907 -2.534215 FALSE FALSE TRUE
## chain_step_id
## 1 1
## 2 2
## 3 3
## 4 4
## 5 5
##
## ...
##
##
## 69.44%
## ====================================================================
## Model 1
## --------------------------------------------------------------------
## Rate parameter period 1 0.0709 (0.0645)
## outdegree (density) 0.0000 (0.0000)
## indegree - popularity 0.0000 (0.0000)
## outdegree - activity 0.0000 (0.0000)
## XW=>X closure of self$component_1_coDyadCovar 0.1000 (0.0000) ***
## self$component_1_coCovar alter 0.0000 (0.0000)
## self$strat_1_coCovar ego 0.5000 (0.0000) ***
## ind. pop.^(1/1) weighted self$strat_2_coCovar 1.0000 (0.0000) ***
## --------------------------------------------------------------------
## Iterations 144
## ====================================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05# saomResDf <- function(res, digits=3)
# {
# df <- data.frame(
# DV=res$effects$name,
# Effect=res$effects$effectName,
# Type=res$effects$type,
# Est=c( NA, round(res$theta, digits = digits)),
# se=c(NA, round(res$se, digits = digits)),
# t=c(NA, round(res$theta/res$se, digits = digits)),
# p=c(NA, round(pt(abs(res$theta/res$se), df=Inf, lower.tail = F) * 2, digits = digits)),
# stringsAsFactors = F
# )
# idx.rate <- grep('rate',df$Effect,T,T)
# df$t[idx.rate] <- NA
# df$p[idx.rate] <- NA
# return(df)
# }
#
# saomResDf(env1$rsiena_model)# ##
# # Create SAOM regression comparison Table
# ##
# saomTable <- function(resList, file=NA, nameMap=NA, digits=3, print=TRUE,
# se.brackets=c('(',')'),p.brackets=c('[',']'),
# lin.comb.test.eval.endow=c('transitive triads'),
# show.sig.stars=TRUE, drop.p.col=TRUE, drop.dv.col=TRUE)
# {
# if (se.brackets[1]=='(') se.brackets[1] <- '\'('
# if ( p.brackets[1]=='(') p.brackets[1] <- '\'('
#
# if (class(resList) != 'list') {
# resList <- list(resList)
# }
# behNames <- c()
# netNames <- c()
# for (res in resList) {
# behNames <- c(behNames, names(res$f$Data1$Behaviors))
# netNames <- c(netNames, names(res$f$Data1$bipartites))
# # netNames <- c(netNames, names(res$f$Data1$nets))
# }
# dvNames <- unique(c(behNames, netNames))
#
# nameList <- list()
# for (dv in dvNames) {
# resDvdf <- data.frame()
# for (res in resList) {
# .effDf <- as.data.frame(res$effects)
# # effTypeNames <- .effDf[grep(dv, res$effects$name,T,T), c('effectName','type')]
# effTypeNames <- .effDf[which(.effDf$name==dv & .effDf$include==T), c('effectName','type') ]
# resDvdf <- unique(rbind(resDvdf, effTypeNames))
# }
# nameList[[dv]] <- unique(rbind(resDvdf, effTypeNames))
# }
#
# dfl <- list()
# for (res in resList) {
# for (dv in dvNames) {
# if (dv %in% res$effects$name) {
#
# # ## results dataframe combining both DVs
# .df <- saomResDf(res, digits=digits)
#
# ##** add linear combination test parameter **
# effdf <- res$requestedEffects
# if( any(str_length(lin.comb.test.eval.endow) > 0) ) {
#
# id.eval <- which( effdf$type == 'eval' & effdf$effectName == lin.comb.test.eval.endow)
# id.endow <- which( effdf$type == 'endow' & effdf$effectName == lin.comb.test.eval.endow)
# ##** Recovering Maintenance Effect 2-sided p-val **
# ##** based on covariance matrix rows for tested parameters **
# #
# lc.est <- .df$Est[id.eval] + .df$Est[id.endow]
# #
# lc.se <- c( sqrt(wt.mat %*% res$covtheta %*% t(wt.mat)) )[1]
# #
# lc.t <- lc.est / lc.se
# lc.p <- pnorm(lc.t, lower.tail = T) * 2 ## two sided
# #
# ##** Wald Test Recovering Maintenance Effect **
# # # ## Alternative Chi-sq (df=1) test to produce p-val
# # # (same as two-sided p-val of norm dist with std.err. above)
# # nparams <- nrow( res$requestedEffects )
# # wt.mat <- matrix(0, 1, nparams)
# # wt.mat[1,id.eval] <- 1
# # wt.mat[1,id.endow] <- 1
# # wt.out <- Wald.RSiena(wt.mat, res)
# # lc.chisq <- wt.out$chisquare
# # lc.df <- wt.out$df
# # lc.1side.z <- c(wt.out$onesided)[1]
# # lc.2side.p <- c(wt.out$pvalue)
# }
# ## Append linear combination test to effects dataframe
# .tmp <- data.frame(
# DV = dv,
# Effect = lin.comb.test.eval.endow,
# Type = 'linear_combination',
# Est = lc.est,
# se = lc.se,
# t = lc.t,
# p = lc.p
# )
# ## append to top if linear combination is in first DV; else append to bottom
# .df <- if (.df$DV[1]==dv){ rbind(.tmp, .df) }else{ rbind(.df, .tmp) }
# ## append to output list
# dfl[[dv]][[ length(dfl[[dv]]) + 1 ]] <- .df[which(.df$DV==dv),]
# }
# }
# }
#
# mod.cols <- c('Est','se','p')
#
# tdf <- data.frame(stringsAsFactors = F)
# for (dv in names(nameList)) {
# hasRowname <- FALSE
# for (rowi in 1:nrow(nameList[[dv]])) { ## effect row
# eff <- nameList[[dv]][rowi,]
# effRow <- list()
# for (modDf in dfl[[dv]] ) { ## model dataframe in DV group
# effId <- which(modDf$Effect == eff$effectName & modDf$Type == eff$type)
# if (length(effId) > 0) {
# effRow[[length(effRow)+1]] <- modDf[effId,mod.cols]
# } else {
# .nadf <- data.frame()
# for (col in mod.cols) .nadf[1,col] <- NA
# effRow[[length(effRow)+1]] <- .nadf
# }
# }
# effRowDf <- cbind(data.frame(DV=dv,Effect=eff$effectName, Type=eff$type, stringsAsFactors = F), cbindDfList(effRow))
# if (!hasRowname) {
# effRowDf <- rbind(effRowDf, effRowDf)
# effRowDf[1,] <- c('', sprintf('Dynamics: %s', dv), rep(NA, ncol(effRowDf)-2))
# hasRowname <- TRUE
# }
# tdf <- rbind(tdf, effRowDf)
# }
# }
#
# # move rate rows to end
# .tmp.rate.row <- tdf[1,]
# .tmp.rate.row$Effect <- 'Rate Parameters'
# rate.idx <- which(tdf$Type=='rate')
# tdf <- rbind(tdf[-rate.idx,], .tmp.rate.row, tdf[rate.idx, ])
#
#
# ##
# obs <- c()
# ns <- c()
# conv <- c()
# convt <- c()
# iter <- c()
# for (res in resList) {
# obs <-c(obs, res$observations)
# ns <- c(ns, attributes(res$f$Data1$nets[[1]][[1]][[1]])$nActors)
# conv <- c(conv, res$tconv.max)
# convt<- c(convt, max(abs(res$tconv)))
# iter <- c(iter, res$n)
# }
#
# # est idx
# idx.est <- which(names(tdf)%in% 'Est')
# #
# tdf[nrow(tdf)+1, ] <- NA
# tdf[nrow(tdf), 'Effect'] <- 'Time Periods'
# tdf[nrow(tdf), idx.est] <- obs
# #
# tdf[nrow(tdf)+1, ] <- NA
# tdf[nrow(tdf), 'Effect'] <- 'Num. Firms'
# tdf[nrow(tdf), idx.est] <- ns
# #
# tdf[nrow(tdf)+1, ] <- NA
# tdf[nrow(tdf), 'Effect'] <- 'Largest converg. t ratio'
# tdf[nrow(tdf), idx.est] <- round(convt, digits = digits)
# #
# tdf[nrow(tdf)+1, ] <- NA
# tdf[nrow(tdf), 'Effect'] <- 'Overall max. converg. ratio'
# tdf[nrow(tdf), idx.est] <- round(conv, digits = digits)
# #
# tdf[nrow(tdf)+1, ] <- NA
# tdf[nrow(tdf), 'Effect'] <- 'Iterations'
# tdf[nrow(tdf), idx.est] <- iter
#
# idx.se <- grep('^se\\.{0,1}',names(tdf),T,T)
# idx.p <- grep('^p\\.{0,}', names(tdf),T,T)
# for (i in 1:length(idx.se)) {
# sei <- idx.se[i]
# pi <- idx.p[i]
# ## Format Standard Error Column
# tdf[,sei] <- apply(tdf[,c(sei,pi)],1,function(x) {
# se <- x[1]
# p <- x[2]
# spfstr <- if(show.sig.stars) {
# sprintf('%s%s%s.%sf%s%s',se.brackets[1],'%',digits+2,digits,se.brackets[2],getPvalStars(p))
# } else {
# sprintf('%s%s%s.%sf%s',se.brackets[1],'%',digits+2,digits,se.brackets[2])
# }
# ifelse(is.na(se)|se=='NA',' ',sprintf(spfstr,as.numeric(se)))
# })
# ## Format p-Value column
# tdf[,pi] <- apply(tdf[,c(sei,pi)],1,function(x) {
# se <- x[1]
# p <- x[2]
# pIsNum <- !is.na(p) & p!='NA'
# spfstr <- sprintf('%s%s%s.%sf%s',p.brackets[1],'%',digits+4,digits,p.brackets[2])
# # cat(sprintf('\np = %s, se = %s, pIsNum = %s, p.spfstr = %s \n',p,se,pIsNum,spfstr))
# ifelse(!pIsNum, ' ',
# ifelse(p<0.001,'[<0.001]',sprintf(spfstr, as.numeric(p))))
# })
# }
#
# ## add Type to name (not eval or rate)
# idx.nonrate <- which(tdf$Type %in% c('endow','creation'))
# tdf$Effect[idx.nonrate] <- apply(tdf[idx.nonrate,c('Effect','Type')],1,function(x){
# sprintf('%s: %s',x[2],x[1])
# })
# tdf <- tdf[,which(names(tdf) != 'Type')]
#
# ## name mapping for effects
# if (!any(is.na(nameMap)))
# {
# ord <- c()
# for (eff in names(nameMap)) {
# idx.nm <- which(tdf$Effect == eff)
# ord <- c(ord, idx.nm)
# if (length(idx.nm) > 0) {
# tdf$Effect[idx.nm] <- nameMap[[eff]]
# }
# }
# tdf <- rbind(tdf[ord,], tdf[-ord,])
# }
#
# if (drop.dv.col) {
# tdf <- tdf[,-1]
# }
# if (drop.p.col) {
# idx.p.col <- grep('^p[\\.\\d]{0,}',names(tdf),T,T)
# tdf <- tdf[, -idx.p.col]
# }
#
# if (print) {
# print(tdf)
# }
#
# if (!is.na(file)) {
# write.csv(tdf, file = file, na = "", row.names = F)
# cat(sprintf('\nFile written to %s\n', file))
# }
#
# return(tdf)
#
# }
# saomTable(list(env1$rsiena_model))# model <- env1$rsiena_model
# # Extract parameter estimates
# estimates <- model$theta
#
# # Extract standard errors from the covariance matrix
# std_errors <- sqrt(diag(model$covtheta))
#
# # Get effect names
# effect_names <- model$effects$effectName
#
# # Compute t-ratios
# t_ratios <- estimates / std_errors
#
# # Compute p-values (assuming normal distribution)
# p_values <- 2 * (1 - pnorm(abs(t_ratios)))
#
# # Function to assign significance stars
# get_significance <- function(p_value) {
# if (p_value < 0.001) return("***")
# else if (p_value < 0.01) return("**")
# else if (p_value < 0.05) return("*")
# else return("")
# }
#
# # Apply significance function
# significance_stars <- sapply(p_values, get_significance)
#
# # Create a results data frame
# results_table <- data.frame(
# Effect = effect_names,
# Estimate = round(estimates, 3),
# `Std. Error` = round(std_errors, 3),
# `t-Ratio` = round(t_ratios, 3),
# Significance = significance_stars
# )
#
# # Print results nicely
# print(results_table)
#
# # For better formatting in RMarkdown:
# library(knitr)
# kable(results_table, caption = "RSiena Model Results", digits = 3)## 1st and last state of the bipartite matrix system
print(env1$bipartite_matrix_init )## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 0 0 0 0 0 0
## [8,] 0 0 0 0 0 0 0 0 0 0 0 0
## [9,] 0 0 0 0 0 0 0 0 0 0 0 0
## [10,] 0 0 0 0 0 0 0 0 0 0 0 0
## [11,] 0 0 0 0 0 0 0 0 0 0 0 0
## [12,] 0 0 0 0 0 0 0 0 0 0 0 0print(env1$bi_env_arr[,, dim(env1$bi_env_arr)[3] ] )## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 0 0 0 0 1 0 1 1 1 1 0 1
## [2,] 0 1 0 1 0 1 0 0 0 1 0 0
## [3,] 0 0 0 0 0 1 0 1 0 1 0 1
## [4,] 1 0 0 0 1 0 1 1 0 0 0 1
## [5,] 0 0 0 0 1 1 1 0 0 0 1 1
## [6,] 0 1 1 1 1 0 0 0 1 1 1 0
## [7,] 0 1 1 0 0 1 1 0 0 0 1 0
## [8,] 1 1 0 1 0 1 0 1 0 0 1 0
## [9,] 1 0 1 0 1 0 1 1 1 1 1 1
## [10,] 1 0 0 1 0 1 1 0 0 1 0 1
## [11,] 0 0 0 0 0 0 1 1 1 1 1 0
## [12,] 0 0 1 1 0 1 0 0 0 0 1 11.4 Print Bipartite Search Network Snapshots
Snapshopts of the the biparite network and the social and component epistasis interactions are taken at fixed step intervals to show the evolving multidimensional coupled search environment (actors and components)
snapshot_ids <- c(1, (1:4)*10, dim(env1$bi_env_arr)[3] )
env1$plot_snapshots(snapshot_ids)
1.5 Plot Actor Degrees (Component Scope and Common Affiliation Social Ties)
Time series is simulated decision steps.
##
env1$plot_actor_degrees(loess_span = 0.5)## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
1.6 Plot Component Degrees (Membership/Popularity and Common Actor Epistasis Ties)
##
env1$plot_component_degrees(loess_span = 0.5)## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
1.7 Plot Utility by Strategy over Time (Simulated Decision Steps)
## [1] 1728 5## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
1.8 Plot Utility by Actor over Time (Simulated Decision Steps)
1.9 Plot Contributions to Actor Utility colored by Strategy
Time series is simulated decision steps.
##
env1$plot_utility_components(loess_span=0.25)## `geom_smooth()` using formula = 'y ~ x'